# Goal: Create Figures 2 & 4 and Table A5
# Jennifer Pan and Yiqing Xu

###########################################
## Stability and Predictability (Figure 2)
###########################################


library(ggplot2)
load("output/s_temporal.RData")
dimnames <-  c("nati", "poli", "econ", "trad", "equi", "ethn")
labels <- c("Nationalism", " Political Liberalism", 
	" Market Economy", " Traditionalism", 
	" Social Equality", " Minority Accommodation")
#width <- rep(4, 6) 
width <- c(9,9,9, 4,4,4)
height <- c(3.5,3.5,3.5,4,4,4)
pos <- c(rep(7, 3), rep(3, 3))
xpos <- c(rep(7.5, 3), rep(4, 3))
ypos <- c(rep(1, 3), rep(1, 3))
size <- c(rep(6, 3), rep(6, 3))
#mycol <- c("gray10","gray50");    
mycol.fill <- c("#E41A1C", "#A6CEE3")
mycol.border <- c("#CB181D", "#045A8D")
labs <- c(" Stability  ", " Predictability  ")
dodge <- 0.95
for (i in 1:6) {
	data <- predict[[i]]
	mat <- data[,c(1,2)]
	colnames(mat) <- c("nitems","corr")
	mat$type <- "a"
	mat.p <- data[,c(1,3)]
	colnames(mat.p) <- c("nitems","corr")
	mat.p$type <- "b"
	mat <- rbind.data.frame(mat, mat.p)
	mat$type <- as.factor(mat$type)
	levels(mat$type) = labs
	mat$corr <- abs(mat$corr)
	nitems <- length(table(mat$nitems))
	p <- ggplot(mat, aes(x=nitems, y = corr)) + 
		geom_violin(aes(fill = type, colour = type), alpha = 0.8, 
			position=position_dodge(dodge), size = 0.5, scale = "width") + 
	  stat_summary(aes(group = type, colour = type), fun.y=median, geom="point", size=2, position = position_dodge(dodge)) +
		theme_bw() + xlab("#Items") + ylab("Correlation") + ylim(0.0,1) + 
		scale_colour_manual(values = mycol.border) + 
		scale_fill_manual(values = mycol.fill) + 
	  annotate(
	    "text",
	    label = labels[i],
	    x = xpos[i],    # horizontal center
	    y = ypos[i],
	    hjust = 0.5,
	    vjust = 1,
	    size = size[i]
	  )
	if (i %in% c(1:3)) {
		p <- p + theme(legend.position = c(0.9,0.15))  		
	} else {
		p <- p + theme(legend.position = c(0.8,0.1)) 
	}
	p <- p + theme(legend.title=element_blank(), legend.background = element_rect(fill="NA")) + 
	guides(fill=guide_legend(keywidth=0.25,keyheight=0.1, default.unit="inch"))
	## save	
	pdf(paste0("graphs/temp_s_",dimnames[i],".pdf"), 
	 	height = height[i], width = width[i])
	print(p)
	graphics.off()	
}


############################################
## heterogeneity on stability (Figure 4)
############################################

rm(list=ls(all=TRUE))
d <- haven::read_dta("data/sample2.dta")
d <- d[which(d$wave2==1),]

getStab <- function(s) {
	out <- rep(NA, 6)
	dims <- c("poli","econ","nati","trad","equi","ethn")
	for (i in 1:6) {
		var <- paste0("index_",dims[i])
		var2 <- paste0("index2_",dims[i])
		out[i] <- cor(s[,var], s[,var2], use = "pairwise.complete.obs")
	}
	return(out)
}

getStab.boot <- function(s, nboot = 1000) {
	n <- nrow(s)
	store <- matrix(NA, 6, nboot)
	for (i in 1:nboot) {
		store[,i] <- getStab(s[sample(1:n, n, replace = TRUE),])
	}
	out <- matrix(NA, 6, 3)
	out[,1] <- getStab(s)
	se <- apply(store, 1, sd)
	out[,2] <- out[,1] - 1.96*se
	out[,3] <- out[,1] + 1.96*se
	rownames(out) <- c("poli","econ","nati","trad","equi","ethn")
	colnames(out) <- c("cor","ci1","ci2")
	return(out)	 
}

stab.hedu <- getStab.boot(d[which(d$eduyr>=16),])
stab.ledu <- getStab.boot(d[which(d$eduyr<16),])
stab.hknow <- getStab.boot(d[which(d$knowledge>=4),])
stab.lknow <- getStab.boot(d[which(d$knowledge<4),])

mycol <- RColorBrewer::brewer.pal(3, "Set1")[c(1,2)]
dim.lab <- c("Political\nLiberalism", "Market\nEconomy","","",
	"Social\nEquality","Minority\nAccommodation")

pdf("./graphs/stability_education.pdf", height = 7, width = 8)
par(mar = c(4,4,1,1))
plot(1, type = "n", xlab = "", ylab = "", 
	axes = FALSE, xlim = c(0.5, 6.5), ylim = c(0,1)); box()
axis(2); axis(1, at = 1:6, labels = FALSE);
axis(1, at = 3, labels = c("Nationalism"));
axis(1, at = 4, labels = c("Traditionalism"));
mtext("Intertemporal Stability", 2, line = 2.5, cex = 1.5) 
axis(1, at = 1:6, tick = FALSE, labels = dim.lab, line = 1, cex.axis = 1)
for (i in 1:6) {lines(c(i,i)-0.1, stab.hedu[i,2:3], lwd = 3, col = mycol[1])}
for (i in 1:6) {lines(c(i,i)+0.1, stab.ledu[i,2:3], lwd = 3, lty = 2, col = mycol[2])}
points(c(1:6)-0.1, stab.hedu[,1], pch = 16, cex = 2, col = mycol[1])
points(c(1:6)+0.1, stab.ledu[,1], pch = 17, cex = 2, col = mycol[2])
legend("bottomleft", pch = c(16, 17), col = mycol, lty = c(1,2),
	legend = c("College or above","Below college"), 
	bty = "n", cex = 1.5, lwd = 3.5)
graphics.off()


pdf("./graphs/stability_knowledge.pdf", height = 7, width = 8)
par(mar = c(4,4,1,1))
plot(1, type = "n", xlab = "", ylab = "", 
	axes = FALSE, xlim = c(0.5, 6.5), ylim = c(0,1)); box()
axis(2); axis(1, at = 1:6, labels = FALSE);
mtext("Intertemporal Stability", 2, line = 2.5, cex = 1.5) 
axis(1, at = 1:6, tick = FALSE, labels = dim.lab, line = 1, cex.axis = 1)
axis(1, at = 3, labels = c("Nationalism"));
axis(1, at = 4, labels = c("Traditionalism"));
for (i in 1:6) {lines(c(i,i)-0.1, stab.hknow[i,2:3], lwd = 3, col = mycol[1])}
for (i in 1:6) {lines(c(i,i)+0.1, stab.lknow[i,2:3], lwd = 3, lty = 2, col = mycol[2])}
points(c(1:6)-0.1, stab.hknow[,1], pch = 16, cex = 2, col = mycol[1])
points(c(1:6)+0.1, stab.lknow[,1], pch = 17, cex = 2, col = mycol[2])
legend("bottomleft", pch = c(16, 17), col = mycol, lty = c(1,2),
	legend = c("High political knowledge","Low political knowledge"), 
	bty = "n", cex = 1.5, lwd =3.5)
graphics.off()

############################################
## Stability by subgroup (Table A5)
############################################


stab <- matrix(NA, 11, 6)
colnames(stab) <- c("poli","econ","nati","trad","equi","ethn")
stab[1,] <- getStab(d)
stab[2,] <- getStab(subset(d, college == 1))
stab[3,] <- getStab(subset(d, college == 0))
stab[4,] <- getStab(subset(d, knowledge >= 4))
stab[5,] <- getStab(subset(d, knowledge <= 3))
stab[6,] <- getStab(subset(d, income >= 6))
stab[7,] <- getStab(subset(d, income <= 5))
stab[8,] <- getStab(subset(d, english >= 4))
stab[9,] <- getStab(subset(d, english < 4))
stab[10,] <- getStab(subset(d, working ==1))
stab[11,] <- getStab(subset(d, working ==0))

# Multiply by 100 to match table scale
stab <- stab * 100

# Add row labels
rownames(stab) <- c(
  "Full sample", 
  "College or above",
  "Below college",
  "High political knowledge",
  "Low political knowledge",
  "Over 8000 RMB",
  "Below 8000 RMB",
  "English proficient",
  "Not proficient",
  "Working or retired",
  "Have not worked"
)

# Convert to data frame for output
stab_df <- as.data.frame(round(stab, 1))

# Reorder columns to match table (Political liberalism → Market economy → Nationalism → Traditionalism → Social equality → Minority accommodation)
stab_df <- stab_df[, c("poli", "econ", "nati", "trad", "equi", "ethn")]

# Write to CSV
write.csv(stab_df, "tables/stability_sub.csv", row.names = TRUE)


